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Abstract: In the field of biomedicine magnetic beads are used for drug delivery and 
to treat hyperthermia. Here we propose to use self-organized bead structures to iso- 
late circulating tumor cells using lab-on-chip technologies. Typically blood flows 
past microposts functionalized with antibodies for circulating tumor cells. Creat- 
ing these microposts with interacting magnetic beads makes it possible to tune the 
geometry in size, position and shape. We develop a simulation tool that combines 
micromagnetics, discrete particle dynamics and fluid dynamics, in order to design 
micropost arrays made of interacting beads. For the simulation of blood flow we use 
the Lattice-Boltzmann method with immersed elastic blood cell models. Paralleliza- 
tion distributes large fluid and particle dynamic simulations over available resources 
to reduce overall calculation time. 
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1 Introduction 

Circulating tumor cells (CTCs) detach from a tumor and can remain in the blood even after the tumor is 
removed. Their presence increases the chance of new tumors developing. It is important to monitor the 
number of CTCs in the blood but their low concentration when compared to normal blood cells makes this 
difficult (one CTC per 5-10 million blood cells). A new and flexible method is required. 

Lab-on-chips with fixed arrays are designed to exploit the properties of the CTCs in order to filter them. 
The most common approaches are mechanical Q based on the size isolation, or using antibodies [4] based 
on the affinity mechanisms. Our proposed chip technology may offer the possibility to switch between affin- 
ity and size capturing through changing the distance between antibody-coated chains using an external field. 
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In this paper we aim at establishing a software environment capable of simulating blood flows inside 
microfluidic devices. We build up a mathematical model describing the blood flow on cellular scale. The 
blood will be considered as a suspension of liquid blood plasma, immersed blood cells and magnetic beads. 
After creating the mathematical model, we implement it into an existing software package ESPResSo (51 . 
We perform a series of simulations focusing on detection of minimal gap sizes for mechanical filters. We 
also present simulations modelling magnetic structures created from magnetic beads under the influence of 
an external magnetic field. 



1.1 Content of the paper 

In Section 2 we present mathematical models describing the fluidic behaviour. First in Section 2.1 we briefly 
explain the dynamics of the fluid. In Section 2.2 we show how general immersed objects are represented. 
The concrete model for blood cell dynamics is explained in Section 2.3 and finally, in Section 2.4 we focus 
on implementation of the cell-cell and cell-boundary interactions. Section 2.5 shows a summary of the 
creation of magnetic chain barriers, which is explained in detail in (T). 
In Section 3 the results of optimal filter gaps are discussed and compared with literature. 

2 Methods 

In this section the description of the blood flow as well as blood cell modeling and magnetic particle dynam- 
ics are covered. We use two different discretizations. First, for the fluid we use the Eulerian fixed grid. On 
this fixed grid, the Lattice-Boltzmann equations are solved, governing the motion of the fluid. Second, for 
the immersed objects we use flexible non-fixed mesh representing the boundary of the object. This method 
is called the Immersed Boundary Method (IBM). The processes on both grids influence each other and they 
are coupled via fluid-object interactions. Combination of these two approaches creates a fluidic simulation 
environment for deformable and rigid objects (i.e. blood cells and magnetic particles, respectively). 

2.1 Lattice-Boltzmann Method (LBM) 

Consider a uniform lattice consisting of square cells, placed over the rectangular three-dimensional domain. 
Instead of solving the Navier-Stokes equations, which solve the conservation equations of macroscopic 
properties, the LBM models the fluid consisting of Active particles. Such particles perform consecutive 
propagation and collision processes over a discrete lattice mesh. The unknown in the LBM is the distribution 
function for Active particles. Macroscopic properties can be recovered by explicit formulas involving the 
unknown distribution function /. The variable /(jc,e,-,f) is the particle density function for the lattice point 
x, discrete velocity vector e, , and time t. We use the D3Q19 (three dimensions with 19 discrete velocities) 
version of the LBM. The governing equations for the LBM, in the presence of external forces, are 

fi (x + 8, , ei , t + 5f ) = fi (x,ei,t)-~ (f (x,e u t)- f- q (x, e,-, t ) ) + F f (x, e t -, f ) 

x 

where 8 t is the time step, x denotes the relaxation time, jf q is the equilibrium function depending on macro- 
scopic variables velocity u(t) and density p(t), and F,- is the external force exerted by the immersed objects 
on the fluid. 



For the description of the fluid dynamics, we use the LBM rather than the discretized Navier-Stokes 
equations. The LBM relies entirely on localized interactions, while the Navier-Stokes equations involve the 
solution of linear systems that couple all points of the grid. The structure of the LBM therefore greatly facil- 
itates parallelizing the calculations. ESPResSo is a publicly available software package that handles particle 
dynamics as well as Lattice-Boltzmann fluid dynamics in a parallel computing environment @. For our de- 
scribed chip technology it is essential to use many CPUs simultaneously in order to simulate a large amount 
of blood cells. Domain decomposition is used to split the simulation environment into cells according to the 
particle position. Parallelization is therefore supported by the fixed grid structure of the Lattice-Boltzmann 
method, the cell system of the ESPResSo package and the storage of data depending on related cell particles. 



2.2 Immersed Boundary Method (IBM) 

Unlike the Lattice-Boltzmann Method the IBM [8] works with non-fixed points. The boundary of each 
suspended object is represented by a set of discrete Lagrangian Immersed Boundary (IB) points that do not 
need to lie on the fluid grid. The IB points move under the influence of forces originated from three (or four) 
sources: 

1. deformation of suspended objects 

2. collisions between objects or between object and wall 

3. coupling between fluid and objects 

4. (magnetic interactions) 



For the motion of IB points we use the following Newton equation 

d% _, 

where m is the mass of the particle, Xj is the position and Fj is the force exerted on the IB point. 

The fluid dynamics equations are discretized on a fixed, uniform Cartesian grid over the entire domain. 
The boundary of each suspended object is represented by a finite set of IB points. Therefore there is no 
need to re-grid even as the RBCs move and deform. This is a significant advantage of the IBM combined 
with the LBM over the traditional methods which need local refinements around the object's boundary and 
thus remeshing every couple of time iterations. For the evolution of the system, the following four steps are 
performed: 



1. The velocities from fluid grid points are transmitted to nearby IB points. 

2. The forces are calculated at each IB point. They result either from elastic behaviour of the objects, or 
from the movement of the object. 

3. The sum of all forces is transmitted from the corresponding IB point to nearby fluid grid points by 
means of a discrete approximate delta-function 



(a) 



(b) 



Figure 1: a) Deformation of RBC in fluid flow b) IB point of RBC interacting with neighbouring fluid grid 
points 

4. The fluid dynamics equations on the grid as well as the motion equations for the suspended objects 
are solved. 

2.3 Cell Dynamics 

2.3.1 Model of a red blood cell (RBC) 

The RBC membrane is a thin shell treated as a two-dimensional sheet. Its elastic properties are described 
in terms of the following elastic moduli [7]: shear modulus K s , area expansion modulus K A , volume expan- 
sion modulus K v and bending modulus K B . When exposed to blood flow, the shape of a red blood cell may 
change dramatically, although both its volume and surface area remain fairly constant due to its lipid-bilayer 
membrane structure. 

To take the mechano-elastic properties of the cells into account, a triangular mesh is created on top of 
the IB points (Fig. [T^). Geometrical entities in this mesh (edges, faces, angles between two faces) are used 
to model the above mentioned elastic moduli Q. 

Volume constraint K v The inner volume of the RBC is held essentially constant at all times. We use a 
very stiff Hookean "spring" approach to describe how a RBC reacts to a change of inner volume. For each 
facet (triangle) we compute a force F v proportional to the change of the volume which is then distributed to 
the three IB points of the facet. 

Area constraint K A Similarly, the surface area of a RBC is constant at all times. Again, we compute and 
distribute a force F A proportional to the change of the area. This force has the effect of (homocentrically) 
shrinking or expanding the face. It is computed locally within each triangle of the mesh (K A ) and globally 
for the whole cell (K A ). 

Shear modulus K s The stretching force F s between two IB points depends on the bond between them. 
Dependence is nonlinear and resembles a neo-Hookean (hyperelastic) behavior. 

Bending force K B We impose a preferred angle between two faces shearing an edge. The departure 
from this angle defines a bending force F B that is distributed to the IB points residing on the particular edge. 
The rest shape resembles a biconcave disc. 



The detailed description of the calibration process is described by Cimrak et al. O. Stretching and 
relaxing tests were performed to validate the cell. Table Q] shows the resulting constants and geometry 
parameters for the red blood cell. 

2.3.2 Model of a circulating tumor cell (CTC) 

CTCs differ in elastic and geometric properties according to the type of cancer. In our case we are con- 
centrating on ephithelial breast cancer cells. Guck et al. (21 use an optical stretcher to extract the optical 
deformability when applying a force on each side of a tumor cell. They create a force between 105 — 333 pN 
on each side with a light power at the front and the back of SOOmW [10]. With a light power of 600mW 
a normal breast epithelial cell (MCF-10) deforms by around 10.5%. A cancerous cell (MCF-7) deforms 
by around 21% and a metastasic cancer cell (modMCF-7) deforms by, at most, 30%. To create our model 
we take the properties of maximum stretching and minimum force (proportional to the lower light power) 
to get the most deformable cell. When this most deformable tumor cell can't pass a filter, then surely also 
less deformable cells will be captured. Table Q] shows the final parameters that fits to the properties of the 
metastasic modMCF-7 cell. 
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Table 1: Parameters of RBC and CTC 



2.3.3 Implementation of blood cells for parallelization 

The simulation of blood cells is problematic with more computer nodes in parallel. The global area con- 
straint Kg and the volume constraint K need information about every IB point of one cell. Particles in 
ESPResSo, in our case mesh nodes, are only stored locally on computer nodes. Therefore the mesh node 
positions of one cell must be transfered to all participating computer nodes to calculate the volume and 
the area of the blood cells before the actual force calculation starts. All other constraints act only between 
nearby mesh nodes and therefore it is no problem in case of parallelization. 
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Figure 2: Validation of metastasic breast cancer cell modMCF-7 with stretching test 



2.4 Cell-Cell/Cell-Boundary Interaction 

Although our simulations run in a rectangular computational box (Fig. |4^), the actual shape of a microchan- 
nel can be formed by arbitrary boundary walls. For the Lattice-Boltzmann fluid the bounce-back rule is 
implemented in ESPResSo. This creates a virtual boundary between two fixed lattice grid nodes. The pro- 
gram can differ between boundary and fluid nodes. At the boundary the direction of velocity is inverted on 
this particular node. 

There is no explicit boundary between the fluid inside and outside of the IB surface. All the interactions 
are communicated via the forces explained in detail in Section [2721 

Collisions between different cells and between cells and walls are handled via particle potentials. In 
ESPResSo we are using a soft-sphere potential with appropriate factors such that no overlapping between 
cells or between cells and walls can occur. If two particles (mesh nodes) of different cells meet they repel 
each other according to the strength of the potential (Fig. 0. This potential is cut with given interaction 
radius (i.e. repelling occurs only when very close together). Although this collision method is computational 
fast it can be problematic in case of e.g. cell aggregation. 



2.5 Magnetic particle dynamics 

The magnetic part of this work has been explained in detail in U. Fig. 0k shows the simulation set-up. 
Viscous blood flows into a pipe filled with soft-magnetic beads. Magnetic charge sheets above and below 
the pipe, with opposite magnetization orientation provide an external magnetic field. 

The gradient force F g (Eqn. Q]) on a bead is given by the negative gradient of the energy of the magnetic 
dipole moment in in the field B. 



F g = V(m-B) (1) 
After applying the magnetic field chains are created in just a few hundred /us under particle interaction 
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Figure 3: Repelling potential cp between different cell nodes as function of position x 



forces Ff. With the moments in of two nearby beads and the distance r we got a formulation (Eqn. [2]) of 
Fj for bead 2 and vice versa for bead 1 because of interacting forces O. In addition friction and collision 
forces are taken into account. 



Different types of CTCs need different distances between the chains for an optimal filtering. Suresh ifTTIl 
shows an overview of different elastic properties of cancer cells. We can use an additional homogeneous 
bias field to tune the arrangement very flexible (also during the filtration process). Particle interaction forces 
causes a larger gap when the total magnetic moment of the beads is higher. Fig. [4}) shows the dependence of 
the gap size on the strenght of the additional field for different particle susceptibilities. With stronger fields 
one can achieve larger gap sizes. 

According to [1] the magnetic forces are strong enough to keep the chain structures static. That is 
why the micromagnetic part, that lasts only a few hundred /us and blood cell dynamics can be computed 
seperately. The simulation run is divided into several parts: 

1. Magnetic particle dynamics to create stable chains 

2. Fixing the chains in space 

3. Fluid dynamics with inserted blood cells 

4. Interaction of blood cells with boundaries (walls (Fig. [5^) and fixed particle chains (Fig. (2t>)) 
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Figure 4: a) Formation of particle chains through the interaction of the drag force and magnetic forces. 
In the picture the magnetic forces are created by the magnetic field gradient F g and magnetic interactions F,. 
The nonuniform magnetic field is created by 2 magnetic charge sheets, b) Variable gap size with different 
susceptibilities and additional external field 

3 Results 

The main idea of this paper was the modeling of blood cells to obtain optimal filter size for mechanical 
capturing. These results should fit to mechanical filters in literature (3] El [321 and further help to improve 
the arrangement of the magnetic bead chains. 

3.1 Minimum gap size 

To obtain the optimal distance between chains for mechanical filtering of circulating tumor cells the elastic 
properties of the blood cells are tested in a special simulation environment (Fig. [5^). Our blood cell models 
must pass a gap, that is variable in size with steps of l/um, with a fluid velocity of 460 /um/s as used at most 
in CTC-chips [4]. Results show that red blood cells can pass the gap even down to 2/um. For our model of 
the metastasic modMCF-7 breast cancer cell the lower limit is 12/um. Seperating the CTCs and the RBCs is 
easily done in the simulations not considered the white blood cells yet, that can have similar diameter than 
the cancer cells. The results also fit to existing mechanical filters like the membrane slot filter of Lu et al. 
0. Their minimum filter gap has a size of 6/jm 

4 Conclusion 

Flexible filters are important to get a high probability of catching cancer cells. Micromagnetic beads are 
established to create chain barriers on demand. Because of the fast creation of stable bead chains we could 
split the magnetic and the cellular calculations. 

To obtain optimal filter gaps we showed how to model blood cells. They are immersed in a Lattice- 
Boltzmann blood flow that discretizes the velocity space. Interaction of the non-fixed mesh grid of the blood 
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Figure 5: a) Simulation of blood flow through variable gap (in this case 6/nm) to obtain minimum gap size 
for each blood cell, b) Circulating tumor cell captured in magnetic bead trap (side and top view) 

cells and the fixed grid lattice is done with the Immersed Boundary Method where the fluid acts on an object 
and vice versa. For the boundary in the fluidic part a bounce-back rule is defined and in the particle part 
potentials help to describe collisions between cells and also walls. 

In this work we have proved that the existing mechanical filters like [3 ] work properly with their minimal 
filter gaps. Further work will include the determination of more specific gaps for different circulating tumor 
cells, just by replacing the elastic properties. Leucocytes need to be investigated too, to get a full simulation 
of blood cell filtration and especially to capture circulating tumor cells. 
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